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Abstract 

The stability of the one-mode nonlinear solutions of the Fermi-Pasta-Ulam - (3 system is 
numerically investigated. No external perturbation is considered for the one-mode exact 
analytical solutions, the only perturbation being that introduced by computational errors 
in numerical integration of motion equations. The threshold energy for the excitation of 
the other normal modes and the dynamics of this excitation are studied as a function of 
the parameter fi characterizing the nonlinearity, the energy density e and the number N 
of particles of the system. The achieved results confirm in part previous results, obtained 
with a linear analysis of the problem of the stability, and clarify the dynamics by which 
the one-mode exchanges energy with the other modes with increasing energy density. In 
a range of energy density near the threshold value and for various values of the number 
of particles N, the nonlinear one-mode exchanges energy with the other linear modes for 
a very short time, immediately recovering all its initial energy. This sort of recurrence 
is very similar to Fermi recurrences, even if in the Fermi recurrences the energy of the 
initially excited mode changes continuously and only periodically recovers its initial value. 
A tentative explanation of this intermittent behaviour, in terms of Floquet's theorem, is 
proposed. 

PACS numbers: 05.45. +b; 63.20.Ry; 63.10. + a 
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1 Introduction 



Since the computer experiment of Fermi- Pasta- Ulam (FPU) [1], many theoretical and numerical 
investigations followed to explain the unexpected results of the experiment [2] - [7]. What one 
expected according to a theorem of Poincare [8] and a theorem proved by Fermi [9] himself in 
his youth, and what was instead observed, has been narrated in several papers (see e.g. [10]) 
in which various aspects of the experiment have been analysed in the framework of the KAM 
theorem [11], ergodic problem [12], statistical mechanics and chaotic behaviour of dynamical 
systems [13, 15]. 

It is well known [16, 17] that, for a periodic FPU-/? chain with an even number of oscillators, 
an initial condition with only a set of excited modes, all having even (resp. odd) indices only, 
cannot lead to the excitation of modes having odd (resp. even) indices. This means that, if 
one considers the set of all modes partitioned in the two subsets of the even and odd modes, 
an initial excitation, completely contained in one of the two subsets, cannot propagate to the 
other. 

There are also other partitions. For example, partitions exist for which a subset contains 
one mode only [18]. More specifically, for each of the modes 

N N N 2 3 
n=— ; — ; — ; -N: -N, (1) 

4 ' 3 ' 2 ' 3 ' 4 ' K J 



(of course when N has the divisibility property required for n in (1) to be an integer) one 
has that, if only one of these mode is initially excited, it remains excited without transferring 
energy to any other mode. 

An important problem is obviously the stability of these one mode solutions (OMS's), since 
it is reasonable to expect some relation between the loss of their stability and the onset of 
chaos in the system. In some sense, the destabilization of the simplest nonlinear modes can 
provide an "upper bound" estimate for the onset of the large scale chaos. The analysis of the 
stability of a generic OMS is very difficult from a mathematical point of view. Only for the 
case n = N/2 there are analytical results ( see for example [19], [20], [18] ) which estimate 
the threshold energy density for the mode to become unstable. In fact the mode n = N/2 is 
the simpler one, since the different components of the perturbation, in modal space, are all 
decoupled and are described by a single Lame equation. 

In this paper we numerically revisit the problem of stability of the OMS's and we present 
the results of a global analysis. We made a numerical study of the stability of the OMS's as 
a function of the number N of particles and of the product e/i, where e is the energy density 
and /i is the parameter of nonlinearity in the Hamiltonian of the system. The analysis is 
based on the numerical integration of the full nonlinear FPU model. No external perturbation 
is considered for the OMS, the only perturbation being that introduced by computational 
errors in numerical integration of motion equations. This simple method works very well: the 
results obtained confirm previous results obtained with a linear analysis of the problem of the 
stability ([19], [20], [18]) and clarify the dynamics with which a nonlinear one- mode exchanges 
energy with the other linear modes, with increasing energy densities. We have found that, in a 
large range of initial excitation energy density, the energy of the OMS's keeps constant for very 
long times, and for short times it is partially transferred to other linear modes; furthermore the 
OMS's corresponding to the two values n = ^ and n = ^N and the OMS's corresponding to 
n = y and n = ^N have the same stability properties. 
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2 The FPU system 



The FPU system is a one- dimensional chain of oscillators (with unit mass), with weakly non- 
linear nearest-neighbour interaction. The nonlinear forces considered are quadratic, cubic and 
broken linear ones. In the present work we consider only the FPU-/3 model (nonlinear cubic 
forces), with periodic boundary conditions. Calling q n and p n the coordinates and the momenta 
of the oscillators, the model is defined by the following Hamiltonian: 

H = H + H 1 , (2) 

being 

Y n y N 

H 0= nE^ + nE(?W" (Ik) 2 , (3) 
A k=l Z k=l 



ffi = lEWi-a) 4 , (4) 
^ k=i 

with q N+1 = q 1 . 

If we introduce the normal coordinates Qk and Pk of the normal modes through the relations: 

N 

Qk = ^2 SkjQj, (5) 
i=i 

N 

Pk = Y. S kjVj, (6) 
1 , . 2nkj Znkj. 

S kj = ^=( B m—+COS—), (7) 

the harmonic energy of mode k is: 

E k = \{P k 2 + ^k 2 Qk 2 ) (8) 
where, in the case of periodic boundary conditions, 

cu 2 fc = 4sin 2 ^. (9) 

We have ojk = &N-k, so that there are only y different frequencies (if we assume N even, for 
simplicity). 

From (3 - 4), the Hamilton equations are obtained in the variables qk and pk, which, in- 
tegrated by standard methods, allow to calculate the normal modes and the energy of each 
mode. 



N 
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3 The one-mode solutions 

Let us first consider the case \x = 0. In this case all normal modes oscillate independently one 
from another and their energies E k are constants of motion. 

In the anharmonic case (// 7^ 0), the normal modes are not independent and the variables Qk 
have not simple sinusoidal oscillations. All the modes are coupled and the differential equation 
for the k-th mode is [18]: 

Qk= F k (Qi, .., Q^), (k — 1, N — 1), (10) 
where F k is the generalized force in normal coordinate space, given by: 



F k (Ql, Qn-l) — —UkQk — ^J^r X] UiUJjUJiC k ijiQiQjQi (11) 

and 

Cijki = —^i+j+k+i + ^i+j-k-i + ^i-j+k-i + Ai_j_ k+ i, (12) 

being A fc = (— l) m for k = mN, if m is a positive integer, and A k = otherwise. 
The nonlinear one-mode solutions correspond to the values of n: 

N N N 2 Ar 3 Ar 
n = — ; — ; — ; -iV; -AT. 13 
4 3 2 3 4 v 7 

From (11) and (12), one deduces [18] that, if only one of these modes is initially excited, it 
remains excited without transferring energy to any other mode. In this case, the equation of 
motion for the excited mode amplitude Q n is: 

zn 2 r\ f^n C nnnn 3 . . 

Qn= Qn ^ ' ^ ^ 

If we assume that at time t — Q„ 7^ and P n = 0, the solution of (14) is: 

Q n (t) — A cn(Q n t, k), (15) 
where Q n and the modulus k of the Jacobi elliptic function cn both depend on A: 



N-l 



n n = cu r ^l + 5 n A 2 , (16) 

(17) 



k 



6 n A* 



\2(1 + s n A*y 



with 5 n = n eC nnnn /2N. 

The solution (15) is periodic with period T n = AK(k)/Vt n where K(k) is the complete elliptic 
integral of the first kind. The energy of the mode is: 

tt. _ l/p2 1 , ,2 n 2 1 ,, U n QnC nnnn 

&n ~ -{^n + U n Q n + /i — ). (18) 

In the next section, the problem of the stability of the OMS corresponding to n = N/2 will be 
analysed. 
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4 The stability of the one- mode solution N/2 



The stability properties of the nonlinear mode N/2 was studied analytically some years ago. 
In [19], the nonlinear mode N/2 was named "the out of phase mode". This expression derives 
from the fact that from (5) and the properties of S k j one has, if the only excited mode is the 
mode N/2: 



qk = 7N { ~ 1)k Qn/2 (19) 



and then: 

q k+1 (t) = -q k (t), k = l,2,...,N. (20) 

In [19], the stability analysis of the out of phase mode (mode N/2) starts from the equations 
of motion for the variables q k . Due to the relations (20), these equations reduce to a single 
equation, describing the anharmonic oscillations of each particle, whose solution is the Jacobi 
elliptic cosine function. Perturbing this solution, and passing to normal modes variables Q k one 
obtains a Lame equation. The stability of the solutions of this equation, which is an example 
of Hill's equation, and then the stability of the mode N/2, is studied with the Floquet theory. 
A numerical analysis shows that the first mode which becomes unbounded, as energy density 
increases, is the mode k = N/2 — 1. A simple approximate formula, valid for large N and fi = 1 
and derived approximating the Hill's matrix with a 3 x 3 matrix, gives for the threshold energy 
density: 

The problem of stability of the mode N/2 was also tackled in [20], in connection with the study 
of tangent bifurcation of band edge plane waves in nonlinear Hamiltonian lattices. In the limit 
of large N, the formula 

«=* = ^-«*» (22) 
* N 3N 2 N 2 1 ; 

is derived for the bifurcation energy density corresponding to the threshold energy density. 
This result is slightly different from the result (21) and the small difference is probably due to 
the rough estimate of the eigenvalue spectrum of the Hill's matrix in [19]. 

The problem of stability of a OMS was reconsidered subsequently in [18] with a detailed 
analysis of the mode N/2. This is the simpler case, because, as we have seen, the different 
components of the perturbation, in modal space, are all decoupled and can be reduced to the 
single Lame equation: 



x r = -uAl+ 12 " A2Cn ^ N/2t] % r, r = l,...,JV-l. (23) 
To obtain in very simple way this equation, we observe that from (2), (3) and (4) one has: 
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Pk= qk+i + qk-i - 2q k + n[(q k +i - qkf - (qk ~ qk-if]- 



(24) 



If the coordinates q k are affected by some error, then the error on the p k , Aft, will be: 

A p k = Ag fe+ i+Ag fc _i-2Ag fc +3/i[(g fe+ i-g fc ) 2 (Ag fc+ i-Ag fc )-(g fe -g fc _i) 2 (Ag fe -Ag fc _i)].(25) 
From (19) and (20) we have: 

4 

(q k+1 - q k ) 2 = (q k - q k _if = — Qn 2 . (26) 
From (25) and (26) we obtain: 

A p k = Aq k+1 + Ag fc _! - 2Aq k + ^Q 2 n [Aq k+1 + Ag fc _x - 2Aq k }. (27) 
Since A q k = Ap k , (27) reads: 

A q k = [1 + ^Q 2 N}[Aq k+1 - 2Aq k + Aq k ^}. (28) 
Passing to modal variables Q k , we finally have: 

A Q k = -u 2 k [l + } AQ k , (29) 

which is equation (23). 

We want to remark, however, that eq. (29) (and then (23)), on which is based the study 
of the stability reported in [18], is the result of a linear analysis and consequently all its 
implications have only local value. 

We outline now the stability analysis of the one-mode N/2, given in [18] . Let be k the 
modulus of the Jacobian elliptic functions. It is related to the product (5 = e/i by the relation: 

Expressing VLn/2 and ij,A 2 /N , as functions of k, and after rescaling time at each fixed k, each 
of the Lame equations (23) can be put in the standard form: 

y" + [a- v{y + 1) k 2 sn 2 (u, k) ] y = 0, (31) 

where the prime superscript denotes differentiation with respect to the new time u = Vt N / 2 t, y 
stands for the generic variable x r and the parameters a and v are defined by the relations: 

a = \(l + Ak 2 )uj 2 r , u(u + \) = \uj 2 . 
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Putting: 

p = sin 2 (7r r/N), (32) 
the two last relations become: 



a = p(l+4k 2 ), u(u + l) = Qp. 



A given pair (p, k) describes different mode numbers r in systems with different N, all having 
the same r/N value, with the same value of the product (3 = ep for the unperturbed solution 
n = N/2. To solve the stability problem for generic values of v, the equation (31) is reduced 
to a Mathieu equation, by approximating the sn 2 by its first-order Fourier expansion. After a 
further time rescaling r = 7r t'/2 K, where K is the complete elliptic integral of the first kind, 
with modulus k, one obtains the canonical form of Mathieu equation: 

+ [a _ 2 gcos(2r)] 2/ = O, (33) 



where: 



a = p(—) 2 (-5 + 4k 2 + 6§), (34) 

71 K 



12 p 
sinh (nK'/K) ' 



g = - -r7T77777Fv ( 35 ) 



K\k) = KVT^k 2 ', (36) 



and E is the complete elliptic integral of the second kind, with modulus k. Both parameters a 
and q depend on k 2 and p. For each fixed p, which corresponds to fix a ratio r/N, changing k 2 , 
namely (5, one can trace a curve in the (q, a) parameter plane. A stability transition happens 
if this curve intersects characteristic curves separating stable from unstable regions for the 
Mathieu equation. 

The main results of the stability analysis reported in [18] are: 

a) for each mode having p > 1/3, there is a threshold value of (3 for instability. Above 
this threshold the nonlinear mode n = N/2 presents an instability causing growth of the mode 
corresponding to p, through parametric resonance; 

b) on the contrary, modes with p < 1/3 (i.e. r/N < 0.196) are always stable in the linear 
approximation, for any energy density of mode n = N/2, so that, perturbations of this mode 
involving only modes with p < 1/3 never lead to instability. These modes, as well as modes 
with p > 1/3, when the product ep is less than the critical value, can grow only if they are 
triggered by the interaction with other modes that are unstable; 

c) for iV > 4 there are always modes with p > 1/3 so that the mode N/2 can never be 
stable for all energy densities. Since the threshold value of (5 is a decreasing function of p, the 
first modes to go unstable, when ep is increased from zero, are the modes r = N/2 — 1 and 
r = N/2 + 1, which have p = cos 2 (tc/N). Therefore, for each (even) number N of particles, 
there is a non-zero value of j3, function of N , below which the nonlinear mode N/2 is stable. 
Since the critical value tends to zero for p — > 1, the critical threshold for the stability approaches 
zero as the number N of particles is increased without limit; 



7 



d) using power series expansion of the various previous formulas, one obtains, for the critical 
value of threshold, the formula: 

A = ^ + 0(iV" 4 ), (37) 

which confirms the N~ 2 dependence found in [19] and in [20] and the numerical values 7r 2 /3 of 
the coefficient of 1/N 2 found in [20]. 



5 Numerical results for the case N/2 

In this section, we present the results of our numerical analysis of the stability of the one-mode 
solution corresponding to n — N/2, as a function of the product e/i and of the number N of 
particles, based on the numerical integration of the full nonlinear FPU model directly in the 
variables q^, Pk- More precisely, we integrate the equations of motion in the variables qk,Pk 
by means of a bilinear symplectic algorithm of the third order, which has been adapted from 
an algorithm employed previously by Casetti [21]. Initial conditions for the variables qk,Pk are 
obtained in the following way. We excite the OMS, at t — 0, always putting Qn/2 7^ and 
Pn/2 = 0. In all the numerical experiments we fix (j, = 0.1 and change the value of the energy 
density e = E N / 2 /N where 

w - l (P 2 ±<? n 2 , ,, u n/2 Qn/2 -, , . 

&N/2 - ^ l/AT/2 + U N/2 V7V/2 + V ' ' ' 

is the energy of the nonlinear one-mode N/2. If we fix the initial value of -Eat/2 (or equiva- 
lently e), the initial value of Qn/2 is obtained from (38) with Pn/2 — 0. Finally, from inverse 
transformations of (5) and (6), the values of qk{0) and Pk{0) are obtained. 

The normal coordinates Qn/2 and Pn/2 of the nonlinear one-mode and the normal coordi- 
nates Qk and Pk of the other normal modes are calculated at fixed time intervals, multiple of 
the integration step. Then, the study of the stability of the OMS is made through the analysis 
of both the time evolution of Qn/2 and Pn/2 and the evolution of the other modes Qk, Pk which 
are generated through computational errors. 

As concerns the numerical integration, we use values of integration time steps At ranging 
from 0.01 to 0.005. For example, the first value is approximately 1/300 of the smallest period 
of oscillation of the atoms in the harmonic case and allows us to obtain a control of the total 
energy E of the lattice, which ensures a relative error AE/E < 10~ 6 . 

To illustrate the various steps of our method of numerical analysis, let us consider the case 
N = 32 and thus the OMS N/2 = 16. In this case C nnnn = 2 and we have the formula (38) for 
the energy of the nonlinear mode. We take a value of the energy density e and integrate the 
equations of motion for the variables qk and pk- The integration time is fixed in such a way to 
observe the instability of nonlinear mode, if the value of (3 — e/j, is greater than the theoretical 
value (37) of the threshold energy density. Typical values of this time are of order 10 6 At where 
At is the integration step. 

For N = 32, we have (3 t = tt 2 /3N 2 = 0.00321. We consider three values of (3: the first one, 
(3 = 0.001, smaller than (3 t) the second one, (3 = 0.005, larger and the third, (3 = 0.1, much 
larger than j3t- In Figs. 1 - 3,the behaviours of the nonlinear mode in the plane Qie,Pm are 
shown for these three values of f3, while in Figs. 4 - 6 we show Q\e(t), as a function of time for 
the same values of (3. 
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Figure 1: P\ 6 vs Qi6 for e/i = 0.001 
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Figure 2: P 16 vs Q w for e/i = 0.005 




Figure 3: P\§ vs Qi6 for e/i = 0.1 
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Figure 4: Q 16 vs t for e/i = 0.001 
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Figure 5: Q\ e vs t for e/i = 0.005 
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Figure 6: Q 16 vs t for e/i = 0.1 
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From the inspection of previous figures, we can observe that, for (3 — e/i very small, well 
below the threshold energy density, the nonlinear OMS N/2 is stable and Qi 6 is a periodic 
function with the same amplitude and the same period of the analytical solution. For e/i = 
0.005, above the threshold, the situation is very different. The period of oscillation is equal 
to the period of the analytical solution, and for very long times, in the plane Qw,Pi6, the 
representative point moves on a one-nonlinear dimensional closed curve, as for values of e/i very 
small; but now, periodically and for short times, the amplitude of oscillation varies, due to a 
decrease of modal energy, and the representative point of the system moves on an open curve 
which tends periodically to shrink. For e/i = 0.1, well above the threshold energy density, we 
observe a chaotic behaviour. This behaviour is also evident from the inspection of the figures 7 
- 9 in which the modal energy is shown, as a function of time, for the three values of (3. In these 
figures, and in all the next figures which show the behaviour of energy vs time, the energy is 
always normalized to the initial excitation energy. 
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Figure 7: Energy of the mode 16 vs t for e/i = 0.001 
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Figure 8: Energy of the mode 16 vs t for e/i = 0.005 



2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 

Timet 



Figure 9: Energy of the mode 16 vs t for eji — 0.1 



As can be seen from the last figures, it is evident that the mode 16 exchanges energy with 
some other mode. According to the theory developed in [19] and in [18], for (3 = 0.005, we are 
above the threshold for the excitation of the adjacent mode 15, which is the first to be excited, 
and the excitation of this mode, due to nonlinear coupling between the modes, should trigger 
other linear modes. In the figures 10 and 11 the behaviour of the the adjacent mode 15 is 
shown. 
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Figure 10: Q 15 vs t for e[i = 0.005 
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Figure 11: Energy of the mode 15 vs t for e/i = 0.005 
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The contemporary excitation of the other modes is shown in Figs. 12 - 14, where the 
energies of the modes n — N/2, n — 15, n — 14 are shown, as functions of time. We remark 
that the energy of the mode 16 is calculated with the formula (18), while, for the other modes, 
the usual formula (8) is utilized. Of course, formula (8) is only indicative for large excitation 
energy, when the variables Qk and lose their meaning of modal variables. 
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Figure 12: Energy of the mode 16 vs t for e/i = 0.005 
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Figure 13: Energy of the mode 15 vs t for e/i = 0.005 
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Figure 14: Energy of the mode 14 vs t for e/i = 0.005 
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The first energy pulse, for the mode 15, appears as in Fig. 15. Near the maximum, the 
value of energy oscillates as in Fig. 16 with a period T ir/2, equal to half the period of the 
oscillations of Q ie . 




Figure 15: Behaviour of the energy of the mode 15 vs t for e/i = 0.005 
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Figure 16: Behaviour of the energy of the mode 15 vs t for e/i = 0.005, near the maximum 
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The initial time interval, necessary to excite, through computational errors, the other modes, 
depends obviously, on the precision of numerical computations. All the previous numerical 
calculations have been performed in double precision. We have observed that, working in 
simple precision, the exchange of energy of the mode N/2 with the other modes occurs much 
before than in double precision. However, since the mechanism is primed, it repeats with the 
same properties either in double or in simple precision. 

The instability of the OMS can also be seen from another point of view, by considering the 
hamiltonian variables We recall that, if the OMS Qn/2 were stable then, from (20), the 
sum qk + qt+i would be always zero. In Fig. 17, as an example, we report the sum g 15 + q ie as 
a function of time. As can be seen from the figure, the instability (sum of the two coordinates 
different from zero) appears when the mode 16 starts to exchange energy with the other modes. 
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Figure 17: q i5 + q ie vs t for e/i = 0.005 when the mode 16 is initially excited. 
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In this context, and to obtain more insight in the behaviour of OMS, it is really very 
interesting to see what are the positions of particles in the chain, i.e. the spatial configuration 
of the chain, in correspondence of a well determined value of energy of the OMS. In figs 18 - 23 
for e/x = 0.005, the values of the coordinates of the 32 atoms, in correspondence of a particular 
value of the energy of mode 16, are shown. To be clear, the representative points of the particles 
are joined by segments. 
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Figure 18: Energy of the mode 16 vs t for e/i = 0.005 




Figure 19: Values of the displacements of atoms in correspondence of the final value of energy 
in Fig. 18 
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Figure 20: Energy of the mode 16 vs t for e/i = 0.005 
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Figure 21: Values of the displacements of atoms in correspondence of the final value of energy 
in Fig. 20 
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Figure 22: Energy of the mode 16 vs t for e/i = 0.005 
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Figure 23: Values of the displacements of atoms in correspondence of the final value of energy 
in Fig. 22 
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As can be seen from the previous figures, the particle chain recovers its " simmetrical form" 
(<Zfc + Qk-i = 0) when the OMS recovers all its initial energy. 



6 The case N/2 as a function of N 

To complete the analysis of the case N/2, in the next figures we show the behaviour of the 
energy of this mode as a function of time for N = 26, 38, 52 and 54, for e/i = 0.005. 
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Figure 24: Energy of the mode N/2 vs t for N = 26 e N = 38. 




Figure 25: Energy of the mode N/2 vs t for N = 52 e N = 54. 



From the previous figures it emerges that, increasing N, the mode N/2 tends to exchange 
all its energy. For N = 52, the mode N/2 periodically loses and recovers roughly all its energy. 
For N > 52, the recovery is not complete and moreover the curve of energy versus time becomes 
irregular. The irregularity and complexity of the curve increase with N and, for N very large, 
the system becomes chaotic. This behaviour is evident in the Figs. 26 and 27, which show the 
behaviour of the energy of the mode N/2, for N = 256 and N = 512. 



21 




8000 



Figure 26: Energy of the mode 128 for N = 256 and e/x = 0.005 
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Figure 27: Energy of the mode 256 for N = 512 and e/x = 0.005 
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Let us now consider the formula (37), obtained analytically, and let us see how we can 
obtain numerically the dependence on N of the threshold energy density. From the analysis of 
the stability reported in [19] and in [18], we know that the first modes to be excited, as the 
energy increases and the mode N/2 becomes instable, are the modes N/2 — 1 and N/2 + 1. 
Then for a value of N, starting from values of (5 — e/i very small, and integrating the motion 
equations for times fairly long, we increase the value of (5 until the first pulse in the energy 
of the mode N/2 — 1 or, equivalently, the first sudden change in the energy of the mode N/2, 
appears. We assume this value of f3 as the value of fit, the threshold energy density. In Fig. 
28, the numerical value of f3 t , determined in this way, and the theoretical value of f3 t (formula 
(37)) are shown for N — 4, 6, . . . , 64. 
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Figure 28: (3 t vs N for the OMS N/2. 
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7 The cases f and §7V 

Following the same procedure used in the case N/2, we have calculated, numerically, the thresh- 
old energy density for the two nonlinear OMS's N/4 and (3/4) N as functions of N. We have 
found that is the same in the two cases. In the case N/4, the initial condition Qn/a 7^ 0, 
Pat/4 = corresponds, in the variables qu, to the initial configuration shown in Fig. 29. 




Figure 29: Initial configuration of the chain particles for the OMS N/4 for N = 32 and 
e/i = 0.01. The displacements from the equilibrium positions are shown vs N. 

We have found that the nonlinear OMS's N/4 and (3/4) JV, for iV = 4, iV = 8 and iV = 12, 
are stable in the limit of an integration time T = 2 x 10 4 . 

The values of (3 t , as a function of N, for the cases N/2 and N/4 are compared in Fig. 30: 
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Figure 30: (3 t vs N for the OMS's N/2 and 7V/4. 



8 The case N/3 and §7V 

In this case, the initial condition: Qn/3 ^ 0, Pn/a = corresponds, in the variables qu, to the 
initial configuration shown in Fig. 31. 
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Figure 31: Initial configuration of the chain particles for the OMS N/3 for N = 30 and 
e/i = 0.008. The displacements from the equilibrium positions, joined by segments, are shown 
vs N. 



The behaviour of (3 t , as a function of N, for the two cases, is given in the Fig. 32: 
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In the following table, for some even value of N between 4 and 64, the threshold (3 t , for each 
OMS, is given. 



N 


(/% 




0% 




(A)|jv 


(A)§jv 


4 


stable 




1.1272 


0.2056 


stable 




6 




stable 


0.1514 


0.0914 




stable 


8 


stable 




0.0665 


0.0514 


stable 




10 






0.0385 


0.0329 






12 


stable 


0.0807 


0.0254 


0.0228 


stable 


0.0807 


16 


0.0447 




0.0136 


0.0128 


0.0447 




18 




0.0250 


0.0106 


0.0101 




0.0250 


20 


0.0219 




0.0086 


0.0082 


0.0219 




24 


0.0138 


0.0128 


0.0059 


0.0057 


0.0138 


0.0128 


28 


0.0096 




0.0043 


0.0042 


0.0096 




30 




0.0079 


0.0037 


0.0037 




0.0079 


32 


0.0072 




0.0033 


0.0032 


0.0072 




36 


0.0056 


0.0055 


0.0026 


0.0025 


0.0056 


0.0055 


40 


0.0045 




0.0021 


0.0021 


0.0045 




42 




0.0040 


0.0019 


0.0019 




0.0040 


44 


0.0037 




0.0018 


0.0017 


0.0037 




48 


0.0031 


0.0031 


0.0015 


0.0014 


0.0031 


0.0031 


52 


0.0027 




0.0013 


0.0012 


0.0027 




54 




0.0025 


0.0012 


0.0011 




0.0025 


56 


0.0023 




0.0011 


0.0010 


0.0023 




60 


0.0021 


0.0021 


0.0010 


0.0009 


0.0021 


0.0021 


64 


0.0018 




0.0009 


0.0008 


0.0018 




66 




0.0018 


0.0009 


0.0008 




0.0018 



The error is of one unit on the last digit. 



In Fig. 33 the behaviour of f3 t as a function of iV is given for the cases N/2, N/3 and N/4. 
The first two values of j3 t of the case N/2 have been left out for clarity. 
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Figure 33: (3 t vs N for the cases N/2, N/3 and N/4. 



9 Discussion 



We have numerically analyzed the stability of the OMS's in the Fermi-Pasta-Ulam system. 
Previously, only for the case N/2 were available theoretical analysis of the stability and ap- 
proximate estimates of the stability threshold for large values of N. Our method is based on the 
numerical integration of full nonlinear differential equations of motion for the particles of the 
chain. The initial conditions for the hamiltonian variables qt and pt are such that only a par- 
ticular nonlinear one-mode is initially excited. No a priori initial perturbation of the analytical 
solution is introduced in the numerical algorithm, the only perturbation being that generated 
by computational errors in numerical integration. Then the method studies the stability of the 
OMS's against the numerical errors introduced by the integration algorithm. We have widely 
analysed the case N/2, which in some sense works as a test, since, for this case, theoretical 
results are available. 

We remark that the OMS's are nonlinear analytical solutions of the complete Fermi-Pasta- 
Ulam system, with linear and nonlinear terms in the hamiltonian, so their stability can't be 
discussed in terms of KAM theorem. The stability of linear modes and nonlinear OMS's, against 
computational errors, are completely disjoint in the sense that a linear mode, initially excited, is 
stable for very long integration times, if the parameter fx is set equal to zero in the hamiltonian. 
Thus the linear modes, excited during the evolution of a nonlinear OMS, are triggered by the 
instability of this nonlinear mode and then only indirectly by the computational errors. This is 
evident in Fig. 34, where the behaviour of energy vs time of the mode N/2 for iV = 32, fx — 0.1 
and e = 0.04 is compared with the behaviour of the energy in the linear case with the same 
value of energy density e = 0.04 and fx — 0. 
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Figure 34: Energy as a function of time for e = 0.04 in the linear case \x = 0, and \x = 0.1. 



As can be seen from this figure, the linear case is stable with respect to the computational 
errors of the numerical algorithm. This different behaviour of the linear mode and of the 
nonlinear OMS is much more evident if we compare, for the same integration time, the orbit in 
the plane Qn/2, Pn/2, for the same value of e = 0.5 and and the two values of /i — 0.0, for the 
linear case, and /i = 0.5 which corresponds, in the nonlinear case, to a value of f3 very much 
larger than the threshold value (3 t = 0.0033. This different behaviour is shown in Fig. 35 and 
in Fig. 36. 
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Figure 35: Orbit in the plane Qi G , P±e for the linear case (ji — 0), with e = 0.5. 



-3-2-10123 

Figure 36: Orbit in the plane Qi6, -Pi6 fc> r the nonlinear case e = 0.5, /i — 0.1. 

Let us first consider the case iV/2. From the inspection of Figs. 1-9, three regimes can 
be observed, varying the product (3 = e\x. For values of (3 very small, well below the threshold 
energy density, the nonlinear OMS N/2 is stable and Qn/2 is a periodic function with the 
same amplitude and the same period of the analytical solution (15). For values of (3 above 
and near the threshold, the situation is very different. The period of oscillation is equal to 
the period of the analytical solution, and for very long times, the representative point moves 
in the plane Qn/2, Pn/2 on a one dimensional closed curve, as for very small values of e/i; but 
now, periodically and for short times, the amplitude of oscillation varies, due to a decrease of 
modal energy, and the representative point of the system moves on an open curve which tends 
periodically to shrink. For e/i = 0.1, well above the threshold energy density, we observe a 
chaotic behaviour. 

In the first regime, where f3 < (3 t , during the whole integration time, the OMS appears 
stable. 
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As soon as f3 exceeds f} t , after an initial time interval, depending on the precision of numerical 
calculations, during which the OMS is stable, the amplification of errors excites the first modes 
which become unstable, namely the modes N/2 — 1 and N/2 + 1. Because of the nonlinear 
coupling between the modes, the excitation of this two modes triggers all the other linear modes. 
The characteristics of this second regime, when the parameter j3 grows, are the increasing 
exchange of energy between the nonlinear mode and the other linear modes and the periodic 
recovery of energy of the nonlinear mode. This regime continues as long as the periodic exchange 
of energy is complete. Further increase of (3, distorts the profile of the curve of energy of the 
nonlinear mode, as a function of time, and subsequently the system becomes chaotic. 

We have seen that, for large values of N, formula (37) gives a good estimate of the instability 
threshold, so we could try to explain the " intermittent behaviour" of the second regime in the 
framework of the theory developed in [18]. As we have seen in section (4), modes with p < 1/3, 
are always stable in the linear approximation, for any energy density of the mode n = N/2, so 
that, perturbation of this mode, involving only modes with p < 1/3, never leads to instability. 
Let us consider for example the case N = 32. We have, from (32), that the modes with r < 7 
are always stable in the linear approximation, for any energy density of mode n = N/2 = 16. 
These modes can grow only if they are triggered by the interaction with other unstable modes. 
As is pointed out in [18], this kind of interaction is neglected in the linear approximation and 
comes into play only when unstable modes have grown and the linearized theory is no longer 
valid. This indirect triggering is evident in the Figs. 37 - 38, where the time evolutions of the 
modes 16 and 6 are compared. No triggering of the mode 6 exists, if (5 < (5 t . 

For (5 > fl t , in the short time intervals, during which the nonlinear one-mode exchanges 
energy with the other modes, the time derivative of the energy of the OMS is not zero and the 
equation (14) and the formula (18) are no longer valid. During these time intervals, a strong 
coupling exists between the mode N/2 and the two adiacent modes and, indirectly, with all the 
other linear modes. However, the exact mechanism by which the nonlinear one-mode loses and 
recovers periodically energy is not clear. We try now to give an explanation of this mechanism 
in terms of Floquet's theorem for parametric oscillators [22], [23]. Let us consider the equation: 

^ + f(t)x = 0, (39) 

where f(t) is a periodic function of time with period T. The Floquet's theorem means that the 
solution of (39) can be written in the form: 

x r (t)=p t J T X r (t) 1 (40) 

where r = 1 or 2, p\p<i — 1 an d X r (t + T) = ±X r (t). The appropriate sign in front of X r (t) is 
determined by the particular form of f(t). With the minus sign one has X r (t + 2T) = X r (t). 
Solution (40) means also that 

x r (t + T) = p r x r {t). (41) 

Thus the values of x r (t), in successive cycles, increases by a factor p when the time interval 
between the observations is equal to the period of f(t). Concerning the stability of the oscillator, 
we can consider three separate cases ([24], [25]). 

The first case corresponds to real values for both pi and pi- Since p\Pi = 1, one of 
them (e.g., p±) is greater than unity. After N — 1 cycles of oscillations, the magnitude of 
the displacement | Xi(0) | grows to the value | p± \ N \ Xi(0) | which exceeds | Xi(0) |: the 
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Figure 37: Energy of the mode 16 vs t for e/i = 0.005. 
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Figure 38: Energy of the mode 6 vs t for e/j, = 0.005. 
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displacement therefore diverges exponentially with iV and the oscillations are unstable. The 
familiar case, known as "parametric resonance", occurs when | /ii |> 1 and the period of the 
Floquet function X(t) is 2T, where T is the period of f{t). 

The second case corresponds to | Hi \—\ 112 \— 1- After N — 1 cycles of oscillation the size 
of the displacement is equal to the initial value, i.e., the oscillations are stable. 

In the third case, n\ — /x| = 1, so that fi± and /i 2 have the same value and the case is 
"degenerate". In this case one can show that the degenerate oscillator is stable or unstable 
with displacement growing linearly with N. 

In the matrix theory of the stability of eq. (39), fi± and ^2 are roots of a quadratic equation. 
If we define the displacement and velocity at the beginning of the n-th cycle of the function 
fit) as x n and v n , respectively, the corresponding variables at the beginning of the next cycle 
of fit) are related to the vector Pin) = (x n ,v n ) by two linear equations which can be written 
as the matrix equation 

P(n + 1) = M(T)P(n). (42) 

Thus the quadratic equation for \i is 

/1 2 - Tr[M(t)}fi + Det[M(T)} = 0, (43) 

where Tr[M{T)} and Det[M{T)] are the trace and determinant of the 2x2 matrix M(T) 
and Det[M(T)\ = = 1. 

The matrix method is particularly useful in determining the stability of the solutions. Sup- 
pose that the graph of f(t) vs t is approximated by a series of steps of width At and that the 
ordinate of the m-th step is 00^. With reference to our case, eq. (23), let us suppose f{t) > 0. 
The transformation matrix [M(m)], relating the displacement x(m + 1) and velocity v (m + 1) 
at the end of the interval At m to the respective values x(m), v (m) at the beginning, is obtained 
by solving the equation of motion, subject to the initial conditions x = x(m) and v = v(m). 
Since f(t) = u)^ n is constant over the interval At m , the equation of motion is that of a simple 
harmonic oscillator so that the following result is easily obtained: 



M(m) = 



I cos(w m At) sin(w m At)/cj r , 
\ — ijj m sin(w m At) cos(cj m At) 



The complete transformation matrix [M(T)], relating the initial and final column vectors, 
is obtained by taking the product of all the matrices for the entire time interval, i.e.: 

M(t) = ]jM(m). (44) 

The question of stability is readily resolved by determining the eigenparameters of M(T), 
that is to say, by solving a quadratic equation. 

The pulsed behaviour of the energy of the nonlinear one-mode solution and the consequent 
pulsed behaviour of the linear modes can then be attributed to the time behaviour of the 
parameters fi r , when f3 > f3 t . Let us consider the Figs. 39 and 40 which show, for N/2 = 16 
and f3 = 0.005, respectively the first energy pulse of the mode 15 and the corrispondent time 
variation of Qi 5 . 
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Figure 39: Energy of the mode 15 vs t for e/i = 0.005. 
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Figure 40: Q 15 vs t for e/i = 0.005. 
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From numerical data we have verified that, in our case, namely in the case of eq. (23), we 
have X r (t + T) = —X r (t) and so X r (t + 2T) = X r (t), where T is the period of the Jacobian 
function cn 2 (Q N / 2 t'i k). For (3 = 0.005, from eqs. (16) and (17), we have Cl N / 2 ~ an d 
cn 2 (QN/2t; k) ~ cos 2 uj^/ 2 t. Since ujn/2 — 2, the period of cos 2 is 7r/2 and so 2T = ir. From 
eq. (41), comparing the values of Qi 5 each interval of 2T = ir, we can obtain the value of 
the parameter /j, r which determines the stability or instability of the solution. The pulsed 
behaviour of the solution and of the energy can then be attributed to a time variation of /i r , 
due to the fact that, during the exchange of energy between the nonlinear one- mode and the 
other linear modes, eq. (23) is no longer exactly valid. If we suppose that the one-mode solution 
Qiq continues to oscillate with its normal frequency, but modulated in amplitude, during the 
energy exchange (see for example Fig. 5), then the parameter \i r varies during this exchange, 
causing the typical pulsed behaviour for (3 > (3 t . Obviously, if (3 is very large, many linear 
modes are involved, the exchange of energy is much greater and irreversible and the OMS does 
not recovers all its initial energy: the whole system tends asintotically to a state of energy 
equipartition. In figure 41 the time behaviour of the parameter fj, is shown for the solution Q i5 , 
when the one-mode excited is the one-mode Qn/2, for N = 32 and (3 = 0.005. With reference 
to Fig. 40, the values of /x are obtained calculating the ratio between two consecutive maxima 
of Qi5- It is clear from Fig. 40 and Fig. 41 the link between the exponential growth and the 
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Figure 41: /ivst for Q l5 , N = 32 and e/i = 0.005. 

subsequent exponential decreasing of Q i5 with the variation of the parameter fj, from values 
greater than one to smaller values. 

10 Conclusions 

In this paper we have studied the problem of stability of the OMS's in the Fermi- Pasta-Ulam f3 
oscillator chain. Although this problem has been tackled many times in the last years, analytical 
results were available only for the case N/2, where N is the number of particles in the chain. 
We have envisaged a simple numerical method, which tests the stability of solutions against the 
numerical errors introduced automatically by the numerical algorithm of integration of motion 
equations. The method reproduces the analytical results already known for case N/2, and 
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allows to obtain the threshold energy density above which the OMS is unstable, in the other 
cases, namely the cases iV/4, N/3, §iV and |iV. 

We have found that, for each case, there is a characteristic energy density, above which, 
there is a large range of values of energy densities, before the chaotic region, in which the 
OMS presents an intermittent behaviour, in the sense that the nonlinear mode keeps its initial 
excitation energy for long times and periodically, abruptly, loses and recovers a fraction of this 
energy. We have verified that, for the case N/2, this characteristic energy density coincides, 
for large N, with the threshold energy density given by (22). Then we have assumed this 
characteristic energy density as threshold energy density also in the other cases. 

We have also obtained that, varying N, the cases y, |iV and ^, |iV have the same 
threshold energy density. 

A tentative explanation of the intermittent behaviour, in terms of Floquet's theorem for 
parametric oscillators, has been given. 
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